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Preface 

My  objective  in  writing  this  thesis  was  to  employ 
fractional  order  time  derivatives  to  build  a  new  model  of 
aeroelastic  behavior.  My  approach  allows  an  intuitive  and 
direct  solution  to  the  equations  of  aeroelastic  instability. 
I  have  restricted  my  development  to  an  airfoil  section  in 
freestream  flow,  so  it  should  provide  a  solid  base  for  more 
complex  models  incorporating  additional  equations  for  three 
dimensional  effects. 

I  wish  to  express  my  heartfelt  thanks  to  LtCol  Ronald 
Bagley.  Without  his  guidance  and  encouragement  I  could  not 
have  begun  this  work. 

I  would  also  like  to  thank  Ted  Fecke  and  the 
propulsion  lab  for  their  support  and  interest  in  this  basic 
development.  A  word  of  thanks  to  the  flight  dynamics  lab  as 
well  for  providing  me  with  some  insight  regarding  current 
applications  of  aeroelastic  modeling. 
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Abstract 

This  thesis  introduces  a  new  model  of  aeroelastic 
behavior.  It  simplifies  the  unsteady  aerodynamic  force 
expressions  by  employing  fractional  order  time  derivatives 
to  create  a  compact  fractional  order  polynomial  that  models 
unsteady  aerodynamic  forces  in  the  transform  domain. 

To  solve  the  flutter  and  divergence  problems,  the 
equations  of  motion  for  a  typical  section  are  expressed  in 
fractional  form.  The  resulting  structural  and  aerodynamic 
equations  of  motion  can  be  combined  into  a  single  fractional 
derivative  eigenvalue  problem.  This  approach  allows  a 
solution  to  the  stability  of  a  lifting  surface  which  obtains 
the  system  eigenvalues  directly  rather  than  relying  on  a 
separately  defined  stability  parameter. 
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The  purpose  of  this  study  was  to  develop  a  closed  form 
approximate  solution  to  the  aeroelastic  instability  problem. 
While  the  intended  application  is  in  turbine  blade  instabil¬ 


ity,  this  study  develops  a  basic  model  not  limited  to 
internal  flow.  The  concepts  presented  are  equally  appli¬ 
cable  to  typical  section  approaches  to  wing  instability  and 
can  be  generalized  to  treat  airfoils  with  control  surfaces. 

A.e.r_oe.  l.a  s  Li  c  i.t^ 

The  typical  section  model  of  an  airfoil  has  a  long 
history  of  use  by  aeroe 1 ast icians .  Most  methods  of 
determining  airfoil  flutter  characteristics  have  their  roots 
in  the  typical  section  model.  Early  work  by  Theodorsen  and 
others  provided  an  exact  solution  for  the  two-dimensional 
flow  of  an  incompressible  fluid  around  an  infinitely  thin 
flat  plate  with  a  wake  composed  of  a  plane  vortex  sheet 
(smooth  trailing  edge  flow)  (17:6-7). 

Theodorsen  went  on  to  develop  a  transcendental  func¬ 
tion  which  describes  the  variation  in  lift  due  to  circu¬ 
lation  resulting  from  the  simple  harmonic  motion  of  an 
oscillating  plate.  Theodorsen’s  function  describes  the 
reduction  and  phase  shift  of  lift  with  increasing  "reduced 
frequency"  (k  =  ub/U  or  actual  frequency  multiplied  by  the 
semi-chord  and  divided  by  the  freestream  velocity).  The 
difficulty  with  this  approach  is  that,  since  Theodorsen’s 


function  is  difficult  to  pose  as  a  time  domain  operator,  it 
has  typically  been  applied  only  to  aerodynamic  forces  resul¬ 
ting  from  pure  sinusoidal  motion.  This  has  helped  various 
"indirect"  methods  to  retain  their  popularity  (6:106). 

Indirect  approaches  to  the  aeroelasticity  problem  rely 
on  the  creation  of  a  "stability  boundary",  with  instability 
presumed  to  occur  when  one  of  the  limits  of  such  a  boundary 
is  exceeded.  While  such  methods  provide  good  estimates  of 
critical  flutter  speeds,  it  is  difficult  to  interpret  their 
results  in  terms  of  stability  margins.  Since  these  methods 
rely  on  the  disappearance  of  a  neutral  stability  solution, 
they  cannot  give  precise  information  about  the  actual 
behavior  of  a  system  away  from  its  critical  points  (3:564- 
568),  (4:235). 

One  such  approach  to  finding  critical  flutter  speeds 
is  the  "U-g"  method.  It  employs  an  artificial  structural 
damping  factor,  g,  which  is  allowed  to  vary  with  airspeed. 
The  "stability  boundary"  for  this  method  is  g  <  0;  when  g  is 
negative,  the  structure  "adds"  energy  to  reach  a  neutrally 
stable  condition.  The  system  becomes  unstable  when  a 
positive  g  would  be  needed  to  produce  a  neutral  stability 
solution  (3:561-568). 

A  more  straightforward  approach  of  solving  for  the 
eigenvalues  of  the  system  has  received  more  attention  in 
recent  years.  The  eigenvalue  solution  provides  a  natural 
stability  limit  at  the  imaginary  axis,  a  simple  solution  to 
the  instability  (flutter  or  divergence)  speeds,  and  a  good 
measure  of  the  degree  of  stability  at  a  given  speed. 
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The  power  previously  demonstrated  by  the  fractional 


calculus  in  modeling  frequency  dependent  properties  and 
describing  linear  phenomena  (2:741-742)  indicated  that  a 
fractional  calculus  based  model  might  produce  an  accurate 
approximation  to  the  reduction  in  lift  and  shift  in  phase 
angle  following  airfoil  movement  modeled  by  Theodorsen’s 
function.  Existing  approximations  to  Theodorsen’s  function, 
such  as  the  Padd  polynomial  and  R.  T.  Jones  approximations 
(8:215,344),  (11:31-38),  (16:6,7)  use  integer  power 

polynomials  which  are  the  Laplace  transforms  of  exponential 
functions  and  their  integer  derivatives.  The  fractional 
derivative,  as  defined  in  equation  (1),  has  a  Laplace  domain 
transformation  which  yields  a  fractional  power  of  the 
Laplace  variable  corresponding  to  the  order  of  the 
derivative  as  shown  in  equation  (2). 


i  d  rt 

Da<  x 1 1 )  >  = 

x(t) 

dT 

(1) 

r(l-a)  dt  J0 

(  t  —  T  )  a 

L{Da<  x ( t ) > )  =  saL{x ( t ) } 

(2) 

A  fractional  calculus  based  "generalized  polynomial" 
would  contain  fractional  powers  of  the  Laplace  variable. 

Such  a  polynomial  would  be  a  Laplace  transform  of  a 
generalized  exponential  function. 

The  model  based  on  this  concept  and  presented  in  this 
work  has  the  advantages  of  accurate  modeling  of  Theodorsen’s 
function  in  a  very  compact  form,  consistent  time-domain 
representations  of  lift,  and  a  direct  eigenvalue  solution. 
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The  first  step  in  development  of  a  fractional 
derivative  model  of  aeroe 1  as t i c i ty  is  the  accurate  modeling 
of  the  Theodorsen  function.  The  Theodorsen’s  function  model 
can  then  be  used  to  develop  the  time  domain  responses  of  the 
system.  The  common  aeroelastic  functions  which  describe 
such  time  responses  are  Wagner’s  and  Kussner’s  functions. 


Theodorsen’s  function  is  a  transcendental  function  of 
reduced  frequency  (k  =  o;b/U)  and  is  an  exact  expression  used 
to  construct  the  circulatory  aerodynamic  forces  resulting 
from  the  harmonic  oscillation  of  a  flat  plate  in  an 
incompressible  flow  (8:212-215).  Various  expressions  are 
available  depending  on  the  desired  argument.  Equation  (3) 
is  an  exact  expression  of  the  Theodorsen  function  on  the 
imaginary  axis.  (This  function  will  operate  throughout  the 
complex  plane  if  ik  is  replaced  by  a  complex  variable,  s.) 

C  ( ik)  =  KiUkl/tKoUM+Kidk)]  (3) 

Kn(ik)  is  a  modified  Bessel  function  of  the  third  kind 
of  order  n  with  an  argument  on  the  imaginary  axis.  The  real 
part  of  the  Theodorsen  function  varies  from  1.0  to  0.5  as 
its  argument  varies  from  zero  to  infinity,  so  modeling  this 
function  using  a  polynomial  requires  that  the  ratio  of  coef¬ 
ficients  of  the  lowest  order  term  must  be  1.0  while  the 
ratio  of  coefficients  of  the  highest  order  term  must  be  0.5. 
Adding  parameters  to  improve  the  fit  of  integer  power 


polynomials  means  adding  higher  powers  of  the  argument.  By 
introducing  a  fractional  power  as  a  parameter,  the  function 
fit  may  be  manipulated  by  reducing  the  order  of  the  poly¬ 
nomial  to  a  fractional  value  as  shown  in  equation  (4). 

C(s)  =  (1  +  Fs0)/(1  +  2Fs0)  (4) 

Optimizations  of  a  fractional  order  polynomial  model 
of  the  Theodorsen  function  over  several  decades  of  reduced 
frequency  result  in  a  fractional  exponent  near  0.82.  Since 
the  optimum  is  not  too  sensitive  to  this  exponent,  a  ratio¬ 
nal  exponent  of  (3  =  5/6  was  chosen  to  simplify  the  treatment 
of  the  resulting  equations.  The  fractional  term  coefficient 
F  was  then  optimized  for  a  wide  frequency  range  resulting  in 
a  value  of  F  =  2.19.  The  fractional  calculus  model  shown  in 
equation  (4)  will  be  shown  to  have  good  accuracy  over  all 
reduced  frequencies  and  a  wide  range  of  complex  arguments. 

Figures  1  through  3  show  the  relative  accuracy  of 
various  approximations  to  the  Theodorsen  function.  The 
accuracy  of  fit  is  compared  to  the  Theodorsen’ s  real  and 
imaginary  parts  over  four  decades  of  reduced  frequency  using 
a  root  mean  square  error  measure. 

Figure  1  shows  the  Theodorsen  function  compared  with 
the  polynomial  model  developed  by  R.T.  Jones  and  used  by  the 
aeroelasticity  module  of  NASTRAN  (16:7)  given  in  equation 
(5).  The  root  mean  square  error  for  this  model  is  0.047. 

N  bn  n  *  0  1  2 

C ( s )  sr  2  bn  =  1  -.165  -.335  (5) 

n=0  l+Bn/s  0n  =  0  .0455  .300 
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Figure  3.  Theodorsen  Function  vs  Fractional  Calculus  Model 


Figure  2  shows  the  Theodorsen  function  vs  a  third 
order  Padd  Polynomial  Model  given  by  equation  (6).  The  root 
mean  square  error  for  this  model  is  0.134. 

(s)3  +  3.5(s)2  +  2.7125(b)  +  .46875 
C(s)  ~  _  (6) 

2 ( s ) 3  +  6 . 5 ( s ) 2  +  4.25(s)  +  .46875 

Figure  3  shows  the  Theodorsen  function  compared  with 
the  fractional  calculus  model  given  by  equation  (4).  The 
root  mean  square  error  for  this  model  is  0.031. 

In  addition  to  providing  an  accurate  model  of  the 
function  for  a  wide  range  of  imaginary  frequencies,  the 
model  should  be  accurate  for  complex  arguments.  Since 
Theodorsen’ s  function  is  composed  of  transcendental 
functions  (Bessel  or  Hankel  functions)  it  has  typically  been 
approximated  by  truncated  series  expansions  or  rational 
polynomials.  Many  of  the  approximations  for  these  functions 
are  useful  only  over  a  limited  range  of  arguments  or  cannot 
be  used  for  complex  arguments. 

Early  work  by  Luke  and  Dengler  (14:478-483)  attempted 
to  expand  Theodorsen’ s  function  to  treat  the  generalized 
oscillations  associated  with  complex  arguments.  While  the 
tables  they  generated  were  submitted  without  proof  and  were 
rejected  in  a  series  of  papers  (18:209),  (13:212),  (12:213) 
much  of  their  work  was  useful.  They  introduced  the  use  of 
R.T.  Jones’  approximation  as  a  model  of  the  generalized 
Theodorsen’ s  function.  Recent  work  by  Edwards  (7:18-24)  has 
provided  a  rigorous  proof  of  the  existence  of  a  generalized 
Theodorsen’ s  function  in  the  form  of  equation  (7)  and  shown 
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that  generalized  aerodynamic  loads  resulting  from  arbitrary 
airfoil  motions  can  be  adequately  treated  using  this 
function. 

C ( s )  =  K1(s)/[K0(S)+Kl(s>]  (7) 

Kn(s)  is  a  modified  Bessel  function  of  the  third  kind 
of  order  n  with  a  complex  argument.  This  indicates  that  an 
approximation  to  Theodorsen’s  function  should  be  able  to 
treat  fully  complex,  rather  than  only  imaginary,  arguments 
over  a  wide  range  of  reduced  frequency.  This  will  allow 
treatment  of  the  stable  behavior  of  an  airfoil  as  well  as 
divergence  and  flutter  predictions. 

Since  the  generalized  Theodorsen  function  operates 
throughout  the  complex  plane,  it  is  more  difficult  to 
approximate  than  the  function  restricted  to  the  imaginary 
axis.  Figure  4  shows  the  generalized  Theodorsen  function 
C(s)  =  F(s)  -  iG(s)  of  argument  I  =  re1®  for  several  values 
of  0  as  calculated  by  Edwards  (7:22).  Figure  5  shows  the 
generalized  Theodorsen  function  calculated  using  the  NASTRAN 
model.  Note  that  the  NASTRAN  approximation  loses  accuracy 
as  it  approaches  the  negative  real  axis  where  its  poles 
occur.  The  imaginary  portion  actually  collapses  back  to 
zero  when  the  argument  lies  on  the  negative  real  axis. 

Figure  6  shows  the  generalized  Theodorsen  function 
calculated  using  the  fractional  calculus  model. 
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Figure  4.  Generalized  Theodorsen  Function 
Edwards  (7:22) 
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For  complex  arguments  near  the  negative  real  axis,  the 
fractional  calculus  model  is  an  excellent  approximation  of 
Theodorsen’s  function.  The  fractional  calculus  model  follows 
the  trends  of  Edwards’  values  throughout  the  complex  plane. 
Even  in  the  areas  approaching  the  negative  real  axis, 
however,  the  fractional  calculus  model  is  more  accurate  than 
other  representations.  Other  models  have  negative  real 
poles  in  the  primary  Laplace  plane  which  cause  more  distor¬ 
tion  near  the  negative  real  axis.  This  overall  accuracy 
allows  the  fractional  calculus  model  to  simultaneously 
predict  the  critical  speed  at  which  an  airfoil  begins  either 
unbounded  oscillations  (flutter)  or  a  monotonic  increase  in 
position  (divergence) . 

Given  that  a  model  of  Theodorsen’s  function  is 
sufficiently  accurate,  it  should  be  possible  to  model  the 
forces  arising  from  an  unsteady  flow.  The  next  step  in 
modeling  Theodorsen’s  function  is  to  compare  the  model  to 
existing  data. 

Experimental  Data 

The  forces  on  an  airfoil  oscillating  harmonically  in  a 
steady  flow  are  described  by  the  Theodorsen  function.  NASA 
Ames  has  published  an  extensive  set  of  experimental  data  for 
oscillating  airfoils  in  subsonic  flow.  The  data  in  this 
section  represent  a  64A010  airfoil  oscillating  in  pitch  in  a 
high  Reynolds  number  subsonic  flow  (5:40). 

The  fundamental  frequency  component  of  the  varying 
dynamic  pressure  at  each  point  along  a  thin  airfoil  can  be 
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interpreted  as  a  complex  number  which  indicates  its  magni¬ 
tude  and  phase  shift  with  respect  to  an  input  motion  (5:11). 
This  complex  pressure  can  be  modeled  using  Theodorsen's 
function.  The  magnitude  of  the  contribution  of  each  point 
to  the  overall  lift  coefficient  varies  along  the  airfoil, 
but  the  pressure  at  each  location  can  be  modeled  using  a 
fractional  calculus  form  of  Theodorsen’s  function  plus  a 
scaled  contribution  to  phase  lag  as  shown  in  equation  (8). 
The  scaling  factors  A  and  d  vary  with  chord  location  (x). 

Cp  =  A(x)C(ik)  +  d(x)(ik)  (8) 

Figure  7  shows  the  variation  in  frequency  dependent 
pressure  coefficient  at  a  typical  chord  location  (x=.033c). 
Figure  8  shows  the  variation  of  the  phase  lag  scaling  fac¬ 
tor,  d,  along  the  airfoil  chord.  Other  pressure  coefficient 
plots  for  other  chord  positions  are  given  in  Appendix  A. 

The  theoretical  expression  for  the  lift  coefficient, 
equation  (9)  (17:12),  depends  on  the  Theodorsen  function  and 
varies  with  reduced  frequency. 

Cj  =  A [C ( ik)  +  0 . 5 ik j  (9) 

This  form  should  result  from  the  integration  of  the 
pressure  coefficients  over  the  entire  airfoil.  Equation  (9) 
indicates  that  the  net  contribution  of  the  phase  lag  scaling 
factor  to  the  overall  lift  coefficient  should  be  0.5. 

Figure  9  shows  the  actual  data  from  the  NASA  test  compared 
to  a  fractional  calculus  model  of  equation  (9). 
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Figure  7.  Complex  Pressure  at  3.3%  chord 
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Figure  8.  Variation  in  Phase  Lag 
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Figure  9.  Lift  Coefficient  Variation 


The  fractional  calculus  model  of  the  Theodorsen 
function  serves  as  a  sound  foundation  for  modeling 
experimental  data.  It  captures  the  theoretical  behavior  of 
the  function  and  so  provides  a  good  model  of  the  physical 
behavior  of  an  airfoil  in  simple  harmonic  unsteady  flow. 
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One  approach  to  applying  Theodorsen’s  function  to  more 
generalized  unsteady  flow  is  to  find  the  time  domain 
response  of  a  typical  section  to  a  prescribed  flow 
condition.  Wagner’s  function  and  Klissner’s  function  are 
common  time  domain  representations  of  the  response  of  a 
typical  section  to  unsteady  flows.  Both  are  dimensionless 
scaling  functions  which  modify  the  lift  which  would  act  on  a 
typical  section  if  the  flow  were  "steady"  at  the  new  values. 
As  a  result  of  this,  both  functions  approach  unity  for  large 
arguments  as  the  lift  builds  to  its  steady  state  value. 

Wagner’s  function  is  the  time  domain  response  of  the 
Theodorsen  function  to  a  unit  step  change  in  circulation. 

It  represents  the  effect  of  a  sudden  change  in  angle  of 
incidence  or  vertical  gust  velocity  on  unsteady  circulatory 
lift.  It  may  be  found  by  solving  the  convolution  of  the 
time  domain  transform  of  Theodorsen’s  function  with  the  unit 
step  function  (3:284). 

The  convolution  which  produces  KUssner’s  function  is 
somewhat  more  complex.  KUssner’s  function  is  the  time 
domain  response  of  Theodorsen’s  function  to  a  sharp-edged 
change  in  vertical  gust  velocity  as  it  is  entered  by  a 
typical  section.  This  results  in  a  convolution  of  the  time 
domain  transform  of  Theodorsen’s  function  with  a  travelling 
wave.  Both  of  these  convolutions  are  performed  in  the 
Laplace  domain  and  transformed  to  provide  time  domain 
expressions . 
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Wagner's  FunctLon 


Wagner’s  function,  k^(a),  describes  the  time  dependent 
circulatory  lift  of  an  airfoil  in  incompressible  flow 
resulting  from  a  step  change  in  circulation.  As  a  result, 
Wagner’s  function  is  a  convolution  of  Theodorsen’s  function 
with  the  unit  step  function.  While  this  approach  is 
conceptually  simple,  it  is  not  possible  to  evaluate  the 
exact  convolution  integral  in  terms  of  simple  functions. 


Wagner’s  function  has  been  approximated  by  convolution 
of  small  argument  expansions  of  the  Bessel  functions  (17:32- 
36) ,  and  has  various  polynomial  and  exponential 
approximations  (8:207).  The  fractional  calculus  model  of 
Theodorsen’s  function,  however,  leads  to  a  simple  and 
accurate  representation  of  Wagner’s  function  through  an 
inverse  Laplace  transform.  In  order  to  find  this 
approximation  of  Wagner’s  function,  it  is  necessary  to 
employ  a  binomial  expansion  on  the  Theodorsen  function's 
fractional  calculus  model.  Equation  (10)  shows  the 
fractional  calculus  model  written  in  a  familiar  form  to 
allow  a  binomial  expansion. 


C(s) 


1  +  FI0 


1  +  2Fse 


(10) 


The  binomial  expansion  leads  to  a  series  in  the  form  of 
equation  (11). 


C(I) 


■  1  -  i  t  1 


n=0 


(ID 
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Each  term  in  the  series  representation  of  C(s)/s  in 


equation  (11)  has  a  known  inverse  Laplace  transform,  so  the 
convolution  integral  may  be  evaluated  by  multiplication  in 
the  Laplace  domain.  Wagner’s  function  can  be  defined  as 
k^(o)  =  L-* {C ( s ) /s ) }  where  a  =  tU/b  is  dimensionless  time  as 
shown  in  equation  (12).  Equation  (13)  gives  the  transformed 
equivalent  in  the  time  domain. 
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The  series  in  equation  (13)  is  a  special  case  of  a 
generalized  exponential  function  called  the  Mittag-Lef f ler 
function  E0(x)  (15:102).  The  definition  of  the  beta  order 
Mittag-Lef f ler  function  is  given  in  equation  (14) 


Efl(x) 


m  r  i 

£ 

n=0  L  T ( 1+nB)  J 


kj^o)  =  u_i(o) 


1  E0 

2 


a 

B-i 

.[.(2F)1/0- 

- 

(14) 


(15) 


Figure  10  shows  the  Mittag-Lef f ler  (order  0=5/6) 
approximation  to  Wagner’s  function  from  equation  (15) 
plotted  against  values  calculated  by  Sears  (17:38)  using 
several  Bessel  function  expansion  equations.  Several 
equations  were  necessary  to  approximate  Wagner’s  function 
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since  tne  expansions  o£  tlie  Bessel  functions  Sears  used 
apply  only  to  limited  ranges  of  the  argument. 


Kiissner’s  function  is  more  complex  than  Wagner’s.  It 
is  a  time  domain  representation  of  the  response  of  an 
airfoil  section  to  a  change  in  circulation  resulting  from  a 
sharp-edged  change  in  vertical  gust  velocity.  It  describes 
the  change  in  circulatory  lift  on  an  airfoil  in  initially 
steady  flow  as  it  penetrates  a  step  change  in  vertical  gust 
velocity  (3:286-289).  As  with  Wagner’s  function,  there  is 
no  way  to  evaluate  the  exact  integral  in  terms  of  simple 
functions.  The  exact  definition  of  Kiissner’s  function  is 
given  in  equation  (16). 

k2(o)  =  r^e-Mdod)  -  I1(s))C(s)  +  I  L  (  s  )  ]  >  (16) 
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By  introducing  the  fractional  calculus  model  of  the 
Theodorsen  function,  this  can  be  rewritten  as  shown  in 
equation  (17).  The  inverse  Laplace  transform  results  in  the 
integral  approximation  given  in  equation  (18). 

k  2 ( o )  =  L-1{e-gl0(s)C(s)  +  e~^I 1 ( s ) [ l-£( s ) ] }  (17) 

fa  d  r-C(g-T)!  T  [u,(t)  -  u, ( t-2 ) ] 

k2(o)  =  _ do _ _ _  dT 

J0  tt[t(2-t)]H 

(l-a)[u_1(o)  -  u_1(o-2)] 

+  _ _ _ _ _  (18) 

u [a ( 2-o) ] H 

Klissner’s  function  can  not  be  readily  simplified  from 
this  form.  Various  other  integral  forms  of  this  function 
are  available  for  numerical  calculations,  but  the  most 
common  representation  used  are  exponential  approximations 
(3:285-288),  (8:344-348). 
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IV.. _ Flutter  Prediction 

In  order  to  predict  system  instability,  we  need  to 
determine  the  system's  eigenvalues.  The  eigenvalues  can  be 
found  by  solving  the  eigenvalue  problem  for  a  typical 
airfoil  section.  Figure  11  shows  a  typical  airfoil  section 
with  its  characteristic  parameters. 

Figure  11.  Typical  Section  Geometry 


For  a  two  degree  of  freedom  model  of  an  airfoil,  the 
Laplace  domain  equations  of  motion  are  written  in  matrix 
form  in  equation  (19).  (7:7-8,165-167) 

{[M]s2  +  [K] }{xs)  =  [L] {xs}/msb3  (19) 


Here 


1  xa 

whJ  0 

r  hs 

CM]  = 

[K]  = 

(xs>  = 

N 

0 

b 

a 

X 

_ j 

_  0  r  a 3  wa 2 

L  a*  j 

The  mass/inertia  matrix  [M]  contains  terms  which 
describe  the  nond imens iona 1 i zed  structural  mass  properties 
of  the  system.  The  lift  and  moment  on  the  section  are 
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coupled  by  xa  which  represents  the  static  imbalance  between 
the  center  of  lift  and  center  of  mass  for  the  section.  The 
ra  term  is  the  nondimens ional i zed  section  radius  of  gyration 
which  represents  the  rotational  inertia  of  the  system. 


The  stiffness  matrix  [K]  contains  terms  which 
represent  the  linear  and  torsional  stiffness  associated  with 
the  airfoil.  These  stiffnesses  are  nondimens ional i zed  to 
produce  the  system’s  natural  frequencies  a)a  and 


The  lift  matrix  [L]  is  more  complex.  It  contains 
terms  which  describe  both  the  circulatory  and  non- 
circulatory  lift  acting  on  an  oscillating  airfoil.  The 
circulatory  terms  are  described  by  Theodorsen’s  function. 
The  nun-c ircu 1 atory  terms  arise  from  the  "apparent  mass" 
lift  which  results  from  the  displacement  of  fluid  by  the 
oscillating  airfoil  (8:446-447). 

[L]  =  {[Mnc]s2  +  { [ Bnc ]  +  C ( s ) [Bc] } (U/b) s 

+  C(s)  [Kc]  (U/b)2}/irM  (20) 


Here 


CM 


ncJ 


-it  ira 

ira  -tt(1/8  +  a) 


[BncJ 


0  -IT 

0  ir(a-J$) 


C  B  c  ] 


2ir 


-1 

a  +  H 


a  -  H 

(a+H) (H-a) 


CKc] 
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(a+H) 


These  equations  can  be  reduced  to  the  eigenvalue 
problem  by  eliminating  the  Laplace  domain  state  vector 
associated  with  the  matrices.  This  results  in  what  would  be 
a  standard  eigenvalue  problem  form  if  the  matrix  [L]  shown 
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in  equation  (21)  were  not  a  function  of  the  Laplace 
variable,  s. 


{[M]s2  +  [K]  -  [L]/msb2){xs}  =  {0}  (21) 

The  circulatory  lift  terms  in  this  equation  are 
described  by  Theodorsen’s  function,  and  the  Bessel  functions 
which  define  Theodorsen’s  function  are  transcendental.  This 
means  that  the  equation  in  exact  form  cannot  be  cast  as  an 
algebraic  eigenvalue  problem  having  constant  matrices 
multiplied  by  integer  powers  of  the  Laplace  variable,  s. 

To  overcome  this  obstacle,  Edwards  employs  an 
iterative  "gradient  search"  method  which  begins  with  an 
initial  guess  for  each  eigenvalue  and  iterates  toward  a 
solution  based  on  a  series  expansion  of  the  Bessel  functions 
in  Theodorsen’s  function.  By  using  the  fractional  calculus 
model  of  Theodorsen’s  function,  constant  coefficient 
matrices  can  be  separated  from  the  frequency  dependent  terms 
in  the  equation.  This  produces  a  set  of  differential 
equations  of  fractional  order.  It  would  be  possible  to  use 
this  approach  with  a  polynomial  approximation,  but  the  poor 
accuracy  of  such  functions  for  complex  arguments  (as  shown 
previously)  would  color  any  results. 
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To  construct  the  fractional  order  eigenvalue  problem, 
it  is  convenient  to  multiply  equation  (21)  by  the 
denominator  of  the  Theodorsen  function  model  as  shown  in 
equat ion  ( 22 ) . 

( l+2Fs5/6) { [M] s2  +  [K]  -  [L] /msb2}{xs}  =  {0}  (22) 

Note  that  the  Theodorsen  function  operates  on  a  non- 
dimensional  Laplace  variable  s  =  sb/U.  It  is  necessary  to 
extract  these  nondimensional izing  factors  from  the 
Theodorsen  function  to  solve  the  characteristic  equation. 
There  is  no  simple  way  to  factor  a  constant  out  of  the 
argument  of  the  exact  representation.  The  fractional 
calculus  model  however,  allows  the  dimensional  terms  to 
appear  explicitly  as  shown  in  equation  (23). 

{ [M] (s2+2Fs17/6(U/b)“5/6) 

+  (l+2Fs5/6(U/b)-5/6){[K]  -  [L]/msb2}{xs}  =  {0}  (23) 

Using  the  definitions  of  these  matrices  given  earlier, 
this  expression  can  be  expanded  to  the  form  given  in 
equation  (24).  (A  more  complete  development  of  this 
equation  is  provided  in  Appendix  B.) 

{s17/6  2F(U/b)-V6{[M]  -  [Mnc]/wp> 

+  s2  {[M]  -  [Mnc]/ir^> 

+  s11/6  (U/b)F{2  [Bnc]  +  [Bc]}/irp 
-  s  (U/b)  { [Bc]  +  [BncJ}/iTM 
+  s5/6  F{2[K] (U/b)~5/6  -  (U/b)7/6  [Kc] /TTp> 

+  [K]  -  (U/b)2[Kc]/*p}{xs}  -  {0} 


(24) 


For  simplification,  the  matrices  associated  with  each 
fractional  power  of  s  will  be  given  the  distinct  notation  of 


equations  (25)  through  (30). 

[ml]  =  2F(U/b)-5/6{[M]  -  [Mnc]/irM>  (25) 

[m2]  =  {[M]  -  [Mnc]/TTM}  (26) 

[m3]  =  ( U/b) F{ 2 [Bnc]  +  [B c]}/up  (27) 

[m4]  =  (U/b)  {  [Bc]  +  [Bnc]}/irM  (28) 

[m5]  =  F{2[K] (U/b)"5/6  -  (U/b)7/6 [Kc] /up}  (29) 

[k]  =  [K]  -  (U/b) 2 [Kc] /up  (30) 

Given  these  matrix  definitions,  equation  (24)  is  rewritten 
{[ml]s^7/6  +  [m2]s2  +  [m3]s^^/6 

+  [m4]s  +  [m5]s5/6  +  [k]}{xs>  =  {0}  (31) 

The  system  eigenvalues  and  the  general  solution  to 


equation  (31)  can  now  be  obtained  directly.  Fractional 
order  differential  equations  can  be  solved  numerically  and 
several  evaluation  techniques,  such  as  the  spectrum  shift 
and  modified  matrix  iteration  techniques,  have  been 
developed  for  numerical  analysis  (9:6-18).  These  methods 
can  be  used  to  solve  the  eigenvalue  problem. 

Rather  than  introduce  these  techniques,  this  thesis 
will  present  an  eigenvalue  solution  using  a  simpler 
algebraic  technique.  This  solution  to  the  fractional  order 
state  equation  has  been  published  (1:488-489)  with  an 
example  for  1/2  order  fractional  powers.  The  penalty  for 
constructing  the  equations  using  this  technique  is  a  higher 
order  equation.  The  solution  for  the  1/6  order  state 
equation  yields  a  pair  of  real ,  symmetric  matrices  of  order 
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34.  These  "pseudomass"  [M]  and  "pseudost i f f ness "  [K] 
matrices  now  form  the  equation  of  motion  for  an  airfoil. 

The  fractional  order  equations  of  motion  are  shown  in 
equation  (32)  where  {X}  is  a  state  vector  containing  all  the 
fractional  derivatives  of  the  heave  and  pitch  coordinates  of 
order  0  to  16/6  in  increments  of  1/6. 


sl/6[Jfl{X}-  +  [K]{X)  =  {0}  (32) 


Here  [M]  = 

0000000000000000  ml 
000000000000000  ml  0 
00000000000000  ml  00 
0000000000000  ml  000 
000000000000  ml  0000 
00000000000  ml  0000  m2 
0000000000  ml  0000  m2  m3 
000000000  ml  0000  m2  m3  0 
00000000  ml  0000  m2  m3  00 
0000000  ml  0000  m2  m3  000 
000000  ml  0000  m2  m3  0000 
0  0  0  0  0  ml  0  0  00  m2  m3  0  0  0  0  m4 
0  0  0  0  ml  0  0  0  0  m2  m3  0  0  0  0  m4  m5 

0  0  0  ml  0  0  0  0  m2  m3  0  0  0  0  m4  m5  0 

0  0  ml  0  0  0  0  m2  m3  0  0  0  0  m4  m5  0  0 

0  ml  0  0  0  0  m2  m3  0  0  0  0  m4  m5  0  0  0 

ml  0  0  0  0  m2  m3  0  0  0  0  m4  m5  0  0  0  0 


[K]  = 

000000000000000  -ml  0 
00000000000000  -ml  0  0 

0000000000000  -ml  0  0  0 

000000000000  -ml  0  0  0  0 

00000000000  -ml  0  0  0  0  0 

0000000000  -ml  0  0  0  0  -m2  0 

000000000  -ml  0  0  0  0  -m2-m3  0 

00000000  -ml  0000  -m2-m3  0  0 

0  0  0  0  0  0  0  -ml  0  0  0  0  -m2-m3  0  0  0 

0  0  0  0  0  0  -ml  0  0  0  0  -m2-m3  0  0  0  0 

00000  -ml  0000  -m2-m3  00000 
0000  -ml  0000  -m2-m3  0000  -m4  0 
000  -ml  0000  -m2-m3  0000  -m4-m5  0 
0  0  -ml  0000  -m2-ra3  0000  -m4-m5  0  0 

0  -ml  0000  -m2-m3  0000  -m4-m5  000 
-ml  0000  -m2-m3  0000  -m4-m5  0000 
0000000000000000k 
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The  eigenvalues  of  the  equation  in  this  form  can  be 
found  using  most  common  matrix  manipulation  programs  without 
using  the  specialized  solvers  mentioned  earlier.  Appendix  C 
contains  a  simple  MATLAB  algorithm  implementation  of  this 
solution.  This  solution  yields  34  eigenvalues  in  the 
Laplace  plane.  Only  two  pairs  of  these  eigenvalues 

(those  on  the  primary  branch)  will  map  into  the  primary 
Laplace  s  plane.  These  eigenvalues  correspond  to  the  pitch 
and  plunge  degrees  of  freedom  of  the  typical  section.  The 
remaining  eigenvalues  fall  on  other  Riemann  sheets  when 
mapped  under  the  z^  transformation  from  the  s^/^  plane  to 
the  s  plane. 

Figure  12  shows  a  typical  distribution  of  eigenvalues 
in  the  s^/6  plane  produced  by  a  MATLAB  matrix  solution  for 
U/b=300.  Each  primary  branch  eigenvalue  in  the  s*/6  plane 
can  be  raised  to  the  sixth  power  to  map  it  onto  the  standard 
Laplace  plane.  Figure  13  shows  the  standard  Laplace  plane 
with  the  primary  branch  eigenvalues  from  Figure  12  mapped 
into  their  proper  locations. 

The  stability  of  the  system  is  easily  determined  since 
the  system  eigenvalues  are  known.  A  complete  root-locus 
plot  of  the  system  eigenvalues  can  be  generated  by  varying 
the  airspeed  factor  U/b.  The  flutter  speed  occurs  when  one 
of  the  system  eigenvalues  crosses  the  imaginary  axis  into 
the  right  half  s-plane.  Divergence  occurs  when  one  of  a 
pair  of  eigenvalues  which  have  reached  the  negative  real 
axis  crosses  onto  the  positive  real  axis. 
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Few  methods  exist  for  determining  aeroelastic  system 
stability  by  actual  extraction  of  system  eigenvalues.  The 
value  of  such  methods  has  been  known  for  some  time.  Goland 
and  Luke  presented  an  early  attempt  at  finding  flutter  speed 
by  iterative  calculation  to  find  aerodynamic  damping 
(10:389-395).  More  recently,  Edwards  has  presented  a  modern 
gradient  search  method  of  iteration  for  section  eigenvalues 
(7:45-51).  In  the  same  paper,  he  introduces  a  Padd 
approximate  and  augmented  states  solution  for  system 
eigenvalues.  This  method  produces  a  state  equation  of 
higher  order  and  requires  considerable  manipulation  of  the 
original  equation  to  reach  a  solvable  form  (7:63-72). 

Figure  14  shows  the  root  locus  plot  generated  by  Ed¬ 
wards  using  his  iterative  solution  (7:51).  Figure  15  shows 
the  root  locus  plot  generated  by  the  fractional  order  state 
equation  solution. 

The  fractional  solution  models  the  behavior  and 
magnitude  of  Edwards’  solution  quite  well.  The  actual 
flutter  speed  of  the  models  varies  slightly  because  Edwards’ 
model  includes  the  effect  of  an  unbalanced  trailing  edge 
flap. 
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Edwards  (7:51) 
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Figure  15.  Root  Locus 


Fractional  Method 


The  fractional  calculus  eigenvalue  problem  formulation 
has  a  number  of  advantages  over  existing  techniques.  It 
provides  a  graphic  and  accurate  analytic  description  of  sys¬ 
tem  stability,  a  measure  of  the  rate  of  approach  to  insta¬ 
bility,  and  a  form  ideal  for  stability  and  control  analysis. 

The  fractional  order  state  equation  method  provides 
several  advantages  over  the  series  expansion  approximations 
of  the  Theodorsen  function.  The  fractional  calculus  model 
is  accurate  over  many  decades  of  reduced  frequency;  the 
series  expansions  require  either  large  numbers  of  terms  or 
more  than  one  series  to  provide  accuracy  for  both  large  and 
small  arguments.  The  fractional  calculus  model  captures  the 
complete  behavior  of  the  function,  but  does  not  produce 
transcendental  characteristic  equations.  This  makes  the 
model  easier  to  manipulate  mathematically  and  allows  the 
nondimensional izing  constants  in  its  argument  to  be  easily 
factored  out.  The  model  has  a  known  inverse  transform,  so 
it  yields  consistent  approximations  to  time  domain 
functions.  These  advantages  over  the  series  approximation 
to  the  Theodorsen  function  allow  the  development  of  an  equa¬ 
tion  of  motion  with  constant  coefficient  matrices. 

The  fractional  order  state  equation  method  also 
provides  several  advantages  over  the  integer  polynomial 
approximations  of  the  Theodorsen  function.  The  polynomial 
models  can  be  used  to  develop  equations  of  motion  with 
constant  coefficient  matrices,  but  problems  arise.  The  poor 
accuracy  of  polynomial  models  for  fully  complex  arguments 
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make  it  difficult  to  apply  them  to  decaying  unsteady 
motions.  The  negative  real  poles  of  integer  polynomial 
models  make  them  difficult  to  apply  to  divergence  instabi¬ 
lity.  These  poles  must  be  shifted  away  from  the  axis  to 
produce  a  divergence  model  in  addition  to  a  flutter  model. 
The  fractional  calculus  model  avoids  this  problem  since  it 
does  not  have  poles  on  the  principal  branch  of  the  s  plane. 

The  fractional  calculus  model  has  fewer  parameters 
than  any  other  available  approximation.  It  reproduces  the 
trends  of  the  Theodorsen  function  with  good  accuracy  for  all 
arguments,  and  so  can  be  employed  in  a  wide  range  of 
unsteady  motions. 

The  advantages  of  the  fractional  order  state  equation 
model  make  it  an  ideal  tool  for  modeling  unsteady  aeroelas- 
ticity.  It  is  a  simple  concept  which  can  be  generalized  to 
more  complex  models.  The  necessary  equations  for  including 
the  effects  of  control  surfaces  are  given  by  Edwards  (7:7— 
8,165-166).  It  could  be  adapted  for  use  in  existing  models 
of  three-dimensional  aeroelastic  behavior,  such  as  turbine 
blade  dynamics  or  helicopter  rotor  dynamics,  which  rely  on 
e i genstructure  determination. 

In  summary,  the  fractional  order  state  equation  model 
has  been  shown  to  be  an  effective  tool  for  aeroelastic  anal¬ 
ysis.  Its  form  is  compatible  with  existing  control  theory 
analysis  techniques  and  will  integrate  easily  with  them.  It 
gives  a  compact,  accurate  and  flexible  approximation  to  the 
effects  of  unsteady  flow  and  allows  direct  solution  of  the 
forces  on  an  airfoil. 
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Apne.il.dix  A:  Complex  Pressure  Coefficient.  Data 


This  appendix  compares  fractional  approximations  of 
unsteady  aerodynamic  pressures  to  data  on  an  airfoil  (NACA 
6 4 A 0 1 0 )  oscillating  in  pitch  about  1/4  chord  at  Mach  0.5  and 
Reynolds  #  =  9 . 9E6  (5:15).  The  real  and  imaginary  parts  of 
the  complex  pressure  are  modeled  using  equation  (8). 

Cp  =  A(x)C(ik)  +  d ( x ) ( i k )  (8) 

Data  Model 


Upper  Surface  Complex 
Real 

A(x=.033)  =  -12 
d  (  x  =  .  0  3  3  )  =  4 

-16 

Re  C  I n 

P 

0  k  0.35 

A ( x= . 09 1 )  =  -8.2 
d(x=.091)  =  0 

-16 

Re  C  I 

P 

4 

0  k  0.35 


Pressure  Coefficient 
Imaginary 


-1 


C 

P 


4  . - .  - . - . 

0  k  0.35 


-1 


C 

P 


4 

0  k  0.35 
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Upper  Surface  Complex  Pressure  Coefficient 


Rea  1 

A ( x  = . 24  3)  =  -5.2 
d(x=. 243)  =  -2.5 

-16 


Re  C 

P 


4 

0  k  0.35 

A ( x  = . 4  0  2 )  =  -3.7 
d(x=. 402)  =-3.5 

- 1  6 


Re  C 

P 


4  . .  -  -  . 

0  k  0.35 

A ( x= . 440 )  =  -2.9 
d(x=. 440)  =  -2.5 

-16 


Re  C 

P 


4  . . . 

0  k  0.35 


Imagi nary 


-1 


Im  C 

P 


4 

0  k  0.35 


-1 


Im  C 

P 


0  k  0.35 


-1 


Im  C  : 
P 


0  k  0.35 
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Upper  Surface  Complex  Pressure  Coefficient 


Real  Imaginary 

A ( x= . 5  84 )  =  -1.7 
d  (  x  =  .  5  8  4  )  =  -3.5 

-i6  -i  ; 

Im  C 

P  ■ 

4  - . . . . '  4  ‘ . ‘  . .  . 

0  k  0.35  0  k  0.35 

A ( x= . 7  3  3  )  =  -1 
d ( x= . 733 )  =  -3 


Re  C 

P 


-16 


Re  C 

P 


Im  C 

P 


4  .  '  4 

0  k  0.35  0  k 

A ( x= . 94  1 )  =  -. 15 
d(x=.941)  =  -1.6 


-16 


Re  C 


-1 


Im  C 


P  : 


0.35  0 


0. 35 


0.35 
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Lower  Surface 
Real 

A(x=.034)  =  13.5 
d ( x= . 034 )  =  -4.5 

16 

Re  C 

P 

-4  '  . 

0  k  0.35 

A(x=.094)  =  8.3 
d ( x= . 094 )  =  -1 

16 

Re  C 

P 


Complex  Pressure  Coefficient 

Imaginary 


1 

Im  C  : 

P 

0  k 

1 

Im  C 

P 


-4  . .  -4 

0  k  0.35  0  k 

A ( x= . 243)  =  5 
d ( x= . 243 )  =  2.5 

16  | . . . . .  1  ; 

Re  C  •  Im  C 

P  '  . . ;  . .  if . 4 .  P 

0  k  0.35  0  k 


0.35 


0.35 


0.35 
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Lower  Surface  Complex  Pressure  Coefficient 


Real  Imaginary  . 

A ( x= .341)  =  3.9 
d ( x= . 34  1 )  =  3.5 

16  1  . . 

Re  C  Im  C 

P  P 

0  k  0.35  0  k  0.35 

A ( x= . 394 )  =  3.2 
d(x=. 394)  =  3.5 

16  1  . 

Re  C  Im  C 

P  P 

-4  .  -4  . . 

0  k  0.35  0  k  0.35 

A( x= . 582  )  =  1.7 
d(x= . 582 )  =  4 

16  . . ;  1  : . . . . . _  ■  p  '-.-"-""'-'  . . 

Re  C  Im  C 

P  P 

0  k  0.35  0  k  0.35 
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Lower  Surface  Complex  Pressure  Coefficient 


Real 

A  (  x  = . 7  3  3 )  =  .9 
d ( x= . 733 )  =  3 

16 

Re  C 

P 

-4 

0  k  0.35 

A(x=.923)  =  .15 
d(x=.923)  =  1.3 

16 

Re  C 

P 


0  k  0.35 


Imaginary 


1 


Im  C 

P 


-4  -  '  . . 

0  k  0.35 


1 


Im  C 

P 


-4 

0  k  0.35 


44 


Given  the  system  eigenvalue  problem  (B.l),  we  must 
first  multiply  all  terms  by  the  denominator  of  the 
Theodorsen  function  model  (l+2Fs^/^).  This  produces 
equation  ( B . 2 ) . 

{[M]s2  +  [K]  -  [L]/(msb2)}{xs}  =  {0}  (B.l) 

{  [M] ( s 2+2Fs 2s 5 /6 ) 

+(l+2Fs5/6){[K]-[L]/(msb2)}  }{xs}  =  {0}  (B.2) 

Since  only  the  terms  in  [L]  contain  C(s)  as  shown  in 
equation  (B.3),  we  can  reform  equation  (B.2)  as  (B.4). 

[L]  =  {[Mnc]s2  +  { [Bnc]  +  C(s)[Bc])(U/b)s 

+  C(s)  [Kc]  (U/b)2}/irp  (B.3) 

{  [M] vs2+2Fs2s5/6)  +  ( l+2Fs5/6) [ K ] 

-  (l  +  2Fs5/6){[Mnc]s2  +  [Bnc]s  }/(ir M)> 

+  (l+Fs5/6){[Bc]s  +  [Kc]/(irM)}  }{xs}  =  {0}  (B.4) 

The  nondimensional izing  constants  in  s  =  sb/U  can  be 
factored  out  to  provide  a  consistent  function  of  a  single 
Laplace  variable  s.  This  results  in  equation  (B.5) 

{  [M] (s2+2Fs2(b/U)5/6s5/6)  +  (l+2F(b/u)5/6s5/6) [K] 

-  (l  +  2F(b/U)5/6s5/6){[Mnc]s2  +  (Bnc]s  }/(irM)} 

+  (l+F(b/U)5/6s5/6){[BcJs  +  [Kc] }  }{xs}  =  {0}  (B.5) 


Carrying  out  the  multiplication  and  grouping  like 
order  terms  reduces  the  equation  to  constant  matrices  and 
fractional  powers  of  s  as  shown  in  equation  (B.6). 

{  s17/6  2F(U/b)"5/6{[M]  -  [Mnc]/irM} 

+  s2  C  [M]  -  [Mnc]/irM> 

+  s11/6(U/b)F{2tBnc]  +  [Bc]}/TrM 
-  s  (U/b)  {  [Bc]  +  [Bnc}>/TTM 
+  s5/6F{2  [K]  (U/b)-5/6  -  (U/b)7/6  [Kc] 

+  [K]  -  (U/b)2[Kc]/TTM  }{xs}  =  (0>  (B.6) 

In  this  form,  the  state  equation  can  be  solved 
directly  using  the  spectrum  shift  or  modified  matrix 
iteration  techniques  or  by  the  matrix  expansion  method 
(1:488-489).  The  "pseudomass"  and  "pseudostiffness" 
matrices  generated  by  the  matrix  expansion  method  are  given 
in  equation  (32).  Appendix  C  is  a  MATLAB  routine  for 
solution  using  this  method. 
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Appendix  C:  MATLAB  Matrix  Construction  and  Solution 


%  Input  routine  for  development  of  Pseudomass  and 
%  Pseudostiffness  matrices 

mu=input (’ Enter  mass  ratio  mu  ’) 
xa=input (’ Enter  x-alpha  ’) 

a=input (’ Enter  cg-rotation  axis  measure  a  * ) 
ra2=input (’ Enter  r  sub  alpha  squared  ’) 
wh=input (’ Enter  omega  sub  h  *) 
wa=input (’ Enter  omega  sub  alpha  ’) 

Ub=input (’ Enter  nondimens i ona 1  velocity  U/b  ’) 


%  Pseudomass  and  pseudostiffness  matrix  generation 


F=2 . 19 ;  %Theodorsen  function  coefficient  ’F’ 

alpha=5/6;  %  C(s)  nondimens iona 1 izing  exponent 


A=[l  +  l/mu  xa-a/mu  ;  xa-a/mu  ra2+(  .  125+a,'2) /mu]  ; 
B= [0  -pi ; 0  pi  * (a-. 5 ) ]  ; 

C=2*pi*[-1  a-.5;a+.5  ( a+ . 5 ) * ( . 5-a) ] ; 

D= [wh”  2  0 ; 0  wa ~  2  *  ra2 ]  ; 

E  =  2  *p i * [0  -1 ; 0  a+ . 5]  ; 


ml=2*F*Ub" (-alpha) *A; 
m2  =A ; 

m3=- (Ub~( 1-alpha)/ (pi *mu ) ) * (2*F*B+F*C) 
m4=-(Ub/(pi*mu) ) * (B+C) ; 
m5=2*F*UbA (-alpha) *D-(F*Ub" ( 2-alpha) / (i 
k=D-Ub''2/(pi*mu)  *E; 


%1 7/6 

term 

%12/6 

term 

( s  *  s 

%1 

.1/6 

term 

% 

6/6 

term 

(s) 

% 

5/6 

term 

% 

0/6 

term 

(1) 

%  Construction  of  Pseudomass  and  Pseudostiffness  matrices 
b 1 =zeros ( 8 ) ; 
bm=zeros (2)  ; 

eml=[bm  bm  bm  ml;bm  bm  ml  bm;bm  ml  bm  bm;ml  bm  bm  bm] ; 

em2=[bm  bm  bm  bm;bm  bm  bm  m2;bm  bm  m2  m3 ; bm  m2  m3  bm] ; 

em3= [m2  m3  bm  bm;m3  bm  bm  bm;bm  bm  bm  bm;bm  bm  bm  m4] ; 

em4=[bm  bm  m4  m5 ; bm  m4  m5  bm;m4  m5  bm  bm;m5  bm  bm  bm] ; 


H=[bl  bl  bl  eml;bl  bl  eml  em2;bl  eml  em2  em3;eml  em2  em3  em4] ; 
bn=zeros (32,2); 


M=[bn  H;ml  bm  bm  bm  bm  m2  m3  bm  bm  bm  bm  m4  m5  bm  bm  bm  bm] ; 
K=[-H  bn ; bm  bm  bm  bm  bm  bm  bm  bm  bm  bm  bm  bm  bm  bm  bm  bm  k] ; 
clear  ml  m2  m3  m4  m5  eml  em2  em3  em4  bl  bm  bn  k  H 


%MATLAB  eigenvalue  solver 
v=eig(-inv(M) *K) 
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